read_srm_export <- function(filename, columns = c("peak_name", "RT.min", "basepeak", "area.cpm", "height.cts", "quantitation")) {
filename %>%
# read excel files
read_excel(sheet = "Integration", skip = 42,
col_names = columns, col_types = rep("text", length(columns))) %>%
as_data_frame() %>%
# remove empty rows
filter(!is.na(peak_name), peak_name != "n.a.") %>%
# convert the relevant numeric columns into numbers
mutate_at(vars(RT.min, area.cpm, height.cts), as.numeric) %>%
# remove useless columns
select(-basepeak, -quantitation) %>%
# add filename info
mutate(file_id = gsub("\\.xls", "", basename(filename))) %>%
select(file_id, everything())
}
# get data
all_data <-
# find all excel files
list.files("Ari_F2_Diatomitedata", recursive = TRUE, full.names = TRUE, pattern = "\\.xls$") %>%
# send them to the read method
lapply(read_srm_export) %>%
# combine the data set
bind_rows() %>%
# pull out sample information
#mutate(sample_id = str_match(all_data$file_id, "TSQ\\d+_GB_(.*)$") %>% { .[,2] }) %>%
# get n replicates
group_by(file_id)
#mutate(n_replicates = length(unique(file_id)))
all_data %>%
ggplot() +
aes(x = peak_name, y = area.cpm, color = file_id) +
geom_point(size = 3) +
theme_bw() +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))
depth_and_rock_info <- read_excel(file.path("metadata", "Diatomite paper OG sample data F2 aromatics 06052018.xlsx")) %>%
filter(!is.na(file_id))
kable(depth_and_rock_info)
| file_id | depth | rock.g | process | TLE.mg | maltenes.mg | intl_std_added.ng | p-Terphenyl-d14_injected.ng |
|---|---|---|---|---|---|---|---|
| TSQ2288_AM_OG126-2_D3666.5 | 3474 | 3.430 | yes | 0.02050 | 0.01432 | 1000 | 7.500 |
| TSQ2166_AM_OG210-2_D3575.5 | 3675 | 6.209 | yes | 0.03000 | 0.02778 | 700 | 5.250 |
| TSQ2191_AM_OG229-2_D3420.5 | 3420 | 5.241 | yes | 0.04718 | 0.03633 | 950 | 7.125 |
| TSQ2163_AM_OG231-2 | 3402 | 5.910 | yes | 0.04357 | 0.02838 | 950 | 7.125 |
| TSQ2291_AM_OG237-2_D3756.5 | 3756 | 6.131 | yes | 0.00287 | NA | 500 | 3.750 |
data_by_depth <-
all_data %>%
left_join(depth_and_rock_info, by = "file_id") %>%
group_by(file_id) %>%
mutate(
n_peaks = n(),
n_standards = sum(peak_name == "d14-pTerph"),
ref_area.cpm = area.cpm[peak_name == "d14-pTerph"],
amount.ug = area.cpm / ref_area.cpm * intl_std_added.ng,
conc_rock.ug_g = amount.ug / rock.g,
conc_TLE.ug_ug = amount.ug/ TLE.mg,
conc_matlenes.ug_ug = amount.ug/ maltenes.mg
#total_area.cpm = sum(area.cpm[peak_name != "d14-pTerph"]),
#area.percent = area.cpm / total_area.cpm * 100
)%>% ungroup() %>%
arrange(file_id, peak_name)
standard <- read_excel(file.path("metadata", "D14 calibration.xlsx")) ###read excel
###calibration curve
standard %>%
ggplot() +
aes(x = Known.ng, y = Measured_area.counts, color = calibration) +
geom_smooth(method = "lm", alpha = 0.5) +
geom_point() +
theme_bw() +
theme(legend.position = "none")
calibrations <-
standard %>%
filter(!is.na(calibration)) %>%
nest(-calibration) %>%
mutate(
fit = map(data, ~summary(lm(`Measured_area.counts`~ `Known.ng`, data = .x))),
coefficients = map(fit, "coefficients"),
intercept = map_dbl(coefficients, `[`, 1, 1),
intercept_se = map_dbl(coefficients, `[`, 1, 2),
slope = map_dbl(coefficients, `[`, 2, 1),
slope_se = map_dbl(coefficients, `[`, 2, 2),
r2 = map_dbl(fit, "r.squared")
)
calibrations %>% select(-data, -fit, -coefficients) %>% knitr::kable(d = 3)
| calibration | intercept | intercept_se | slope | slope_se | r2 |
|---|---|---|---|---|---|
| sept2017 | -131054.9 | 2131453 | 2275077 | 90941.37 | 0.997 |
| oct2017 | -1488460.7 | 2799385 | 1557777 | 107323.96 | 0.995 |
These numbers are not useful for anything else.
calib_data <-
data_by_depth %>%
# temp
mutate(calibration = "oct2017") %>%
left_join(calibrations, by = "calibration") %>%
mutate(
total_volume.uL = 100,
total_inject.uL = 1.5,
ref_amount_inject_expected.ng = intl_std_added.ng/total_volume.uL * total_inject.uL * 1000,
ref_amount_inject_measured.ng = (ref_area.cpm - intercept)/slope,
ref_amount_measured.ug = total_volume.uL/total_inject.uL * ref_amount_inject_measured.ng * 1/1000,
yield = ref_amount_inject_measured.ng/ref_amount_inject_expected.ng
)
#calib_data_corr <- calib_data %>% mutate(conc_rock.ug_g = amount_yield_corr.ug / `rock `) #in microgram lipid per gram rock
#View(calib_data_corr)
calib_data %>%
select(file_id, yield) %>%
arrange(file_id) %>%
unique() %>%
ggplot() + aes(file_id, y = 100*yield) +
geom_point(size = 3) +
theme_bw() + theme(axis.text.x = element_text(angle = 90, hjust = 0, vjust = 0.5))
# functions to make it easy to sum up peaks
sum_peaks <- function(df, filter_condition, new_peak_name) {
filter_condition <- sprintf("(%s)", str_c(filter_condition, collapse = "|"))
filter(df, str_detect(peak_name, filter_condition)) %>%
summarize(
file_id = file_id[1],
depth = depth[1],
conc_rock.ug_g = sum(conc_rock.ug_g)
) %>%
mutate(peak_name = new_peak_name)
}
ratio_peaks <- function(df, filter_top, filter_bottom, new_peak_name) {
filter_top <- sprintf("(%s)", str_c(filter_top, collapse = "|"))
filter_bottom <- sprintf("(%s)", str_c(filter_bottom, collapse = "|"))
filter(df, str_detect(peak_name, filter_top) | str_detect(peak_name, filter_bottom)) %>%
summarize(
file_id = file_id[1],
depth = depth[1],
ratio = sum(conc_rock.ug_g[str_detect(peak_name, filter_top)]) / sum(conc_rock.ug_g[str_detect(peak_name, filter_bottom)])
) %>%
mutate(peak_name = new_peak_name)
}
final_data <- data_by_depth %>%
group_by(file_id) %>%
do({
bind_rows(.,
sum_peaks(., "C15", "all C15"),
sum_peaks(., "C16", "all C16"),
sum_peaks(., "C17", "all C17"),
sum_peaks(., "C18", "all C18"),
sum_peaks(., "C19", "all C19"),
sum_peaks(., "C20", "all C20"),
sum_peaks(., "C21", "all C21"),
sum_peaks(., "C24", "all C24"),
sum_peaks(., "C26", "all C26"),
sum_peaks(., "Aryl", "all_Aryl_Isop"),
#sum_peaks(., "MP", "3ring_MP"),
sum_peaks(., c("Acenapthene", "Flourene"), "2ring_all" ),
#sum_peaks(., "Acenapthene", "2ring_1"),
#sum_peaks(., "Flourene", "2ring_2"),
sum_peaks(., c("Phenanthrene", "Flourantherene", "Retene", "MP"), "3ring_all" ),
#sum_peaks(., "Phenanthrene", "3ring_1"),
#sum_peaks(., "Flourantherene", "3ring_2"),
#sum_peaks(., "Retene", "3ring_3"),
sum_peaks(., c("Pyrene", "Benzo[a]anthracene", "Triphenylene", "Chrysene", "flouranthrene"), "4ring_all" ),
#sum_peaks(., "Pyrene", "4ring_1"),
#sum_peaks(., "Benzo[a]anthracene", "4ring_2"),
#sum_peaks(., "Triphenylene", "4ring_3"),
#sum_peaks(., "Chrysene", "4ring_4"),
#sum_peaks(., "flouranthrene", "4ring_5s"),
sum_peaks(., c("Benzo[e]pyrene", "Benzo[a]pyrene", "Perylene", "Ideno", "Dibenz"), "5ring_all" ),
#sum_peaks(., "Benzo[e]pyrene", "5ring_1"),
#sum_peaks(., "Benzo[a]pyrene", "5ring_2"),
#sum_peaks(., "Perylene", "5ring_3"),
#sum_peaks(., "Ideno[c,e]", "5ring_4"),
#sum_peaks(., "Dibenz[a,h]", "5ring_5"),
sum_peaks(., c("Benzo[ghi]", "Coronene"), "6ring_all" ),
sum_peaks(., c("Acenapthene", "Flourene", "Phenanthrene", "Flourantherene", "Retene", "Pyrene", "Benzo[a]anthracene", "Triphenylene", "Chrysene", "flouranthrene", "Benzo[e]pyrene", "Benzo[a]pyrene", "Perylene", "Ideno[c,e]", "Dibenz[a,h]", "Benzo[ghi]", "Coronene"), "all_PAH")
#sum_peaks(., "Benzo[ghi]", "6ring_1"),
#sum_peaks(., "Coronene", "6ring_2")
) }) %>% ungroup()
subset(final_data, peak_name=='d14-pTerph') %>%
ggplot() +
aes(x = amount.ug, y = depth) +
geom_point() +
#facet_wrap(~peak_name, scales = "free") +
scale_y_reverse()
Iso <- subset(final_data, peak_name=='Isorenieretane') %>%
ggplot() +
aes(x = depth, y = conc_rock.ug_g) +
geom_point() +
geom_line() +
#facet_wrap(~peak_name, scales = "free")
scale_x_reverse() +
coord_flip()
ggplotly(Iso)
## We recommend that you use the dev version of ggplot2 with `ggplotly()`
## Install it with: `devtools::install_github('hadley/ggplot2')`
Chlor <- subset(final_data, peak_name=='Chlorobactane') %>%
ggplot() +
aes(x = depth, y = conc_rock.ug_g) +
geom_point() +
geom_line() +
#facet_wrap(~peak_name, scales = "free")
scale_x_reverse() +
coord_flip()
ggplotly(Chlor)
## We recommend that you use the dev version of ggplot2 with `ggplotly()`
## Install it with: `devtools::install_github('hadley/ggplot2')`
subset(final_data, peak_name %in% c('Chlorobactane', 'Isorenieretane')) %>%
ggplot() +
aes(x = depth, y = conc_rock.ug_g) +
geom_point() +
geom_line() +
facet_wrap(~peak_name, scales = "free") +
scale_x_reverse() +
coord_flip()
subset(final_data, peak_name=='all_Aryl_Isop') %>%
ggplot() +
aes(x = depth, y = conc_rock.ug_g) +
geom_point() +
geom_line() +
facet_wrap(~peak_name, scales = "free") +
scale_x_reverse() +
coord_flip()
subset(final_data, peak_name %in% c('all_PAH', '2ring_all', '3ring_all', '4ring_all', '5ring_all', '6ring_all')) %>%
ggplot() +
aes(x = depth, y = conc_rock.ug_g) +
geom_point() +
geom_line() +
facet_grid(~peak_name, scales = "free") +
scale_x_reverse() +
coord_flip()